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Heat transport at nanoscales in semiconductors is investigated with a statistical method. The 
Boltzmann Transport Equation (BTE) which characterize phonons motion and interaction within 
the crystal lattice has been simulated with a Monte Carlo technique. Our model takes into account 
media frequency properties through the dispersion curves for longitudinal and transverse acoustic 
branches. The BTE collisional term involving phonons scattering processes is simulated with the 
Relaxation Times Approximation theory. A new distribution function accounting for the collisional 
processes has been developed in order to respect energy conservation during phonons scattering 
events. This non deterministic approach provides satisfactory results in what concerns phonons 
transport in both ballistic and diffusion regimes. The simulation code has been tested with silicon 
and germanium thin films; temperature propagation within samples is presented and compared 
to analytical solutions (in the diffusion regime). The two materials bulk thermal conductivity is 
retrieved for temperature ranging between 100 K and 500 K. Heat transfer within a plane wall with 
a large thermal gradient (250 K-500 K) is proposed in order to expose the model ability to simulate 
conductivity thermal dependence on heat exchange at nanoscales. Finally, size effects and validity 
of heat conduction law are investigated for several slab thicknesses. 

PACS numbers: 05.10.Ln, 44.10.+i, 63.20.-e, 65.40.-b 



I. INTRODUCTION 

The development of nanotechnologies has lead to an 
unprecedented size reduction of the electronic and me- 
chanical devices. For example, transistors of a few 
nanometer size are now openly considered^. The heat 
that will be dissipated by joule effect in these semi- 
conductors junctions will reach soon the levels of the heat 
dissipated in a light bulb. This high volumetric heat dis- 
sipation in electronic devices will have to be evacuated 
very efficiently in order to avoid possible failures of the 
systems. This task will not be achieve without a sharp 
knowledge of the phenomena governing the heat transfer 
at nanoscale. Furthermore, new technologies based on a 
local heating are being developed in order to enlarge the 
computer hard disk capacity. The ultimate limit of stor- 
age is to write a byte at the atomic scale. This goal is 
already feasible with near-field microscope probes but at 
a too slow rate. A way to write bytes at the nanometer 
scale is the melting of a polymer by heating it on a very 
short time scale (< 1 ns) by an array of heated near-field 
probes^. In this example, the heat transfer has to be con- 
trolled not only at the nanometer scale but also at the 
nanosecond scale. 

Through these two examples, one can anticipate that 
the foreseeing technological challenges in miniaturization 
will have to solve more and more problems of heat trans- 
fer at short time and space scale. However, the physics 
of heat transfer usually used (Fourier's law, Radiative 
Transfer Equation) can no longer be applied when some 
characteristic length scales are reached^. In thermal radi- 
ation for example, wave effects appear as the system char- 
acteristic lengths becomes lower than the typical wave- 
length (A ~ 10^m at T = 300 K)i£. A substantial in- 



crease of the radiative heat transfer can even be reach 
at nanometric distances^. On its side, heat conduction 
is classically described by the Fourier law and the heat 
conduction equation (dT/dt = aAT) which is a diffusion 
equation. It is well known that this kind of equation can 
be interpreted as a random walk of particles^. In the 
case of heat conduction, we are actually dealing with en- 
ergy carriers which are electrons in metals and phonons 
in crystalline materials. When these carriers undergo a 
large number of collisions, the use of the diffusion equa- 
tion is valid whereas a more careful study is required 
when the number of interactions between carriers lowers. 

A way to achieve this goal is to consider the evolu- 
tion of a distribution function /(r,p,i) which describes 
the number of particles in a certain elementary volume 
d 3 rd 3 p around the point (r, p) in the phase space. The 
evolution equation of /, called the Boltzmann Transport 
Equation (BTE) makes / to vary in space and time under 
the influence of advection, external force and collision 8 . 
Note that this approach is not relevant to treat the wave 
aspects of the problem such as interference or tunneling. 
The understanding and the modeling of the collision term 
is actually the key point in the resolution of the BTE. It 
can sometimes be fully modelized as in radiation trans- 
fer. Then the collision term is in that case the sum of an 
absorption term, an elastic scattering term and an emis- 
sion term proportional to an equilibrium distribution 9 . 
Many resolution technique have been developed in radi- 
ation transfer such as the Discrete Ordinates Method, 
the Monte Carlo method or the ray-tracing method 10 . 
They can hardly be used when the collision processes are 
inelastic as it is the case for electrons and phonons 11 . 
For example, the phonons, which are eigenmodes of the 
harmonic oscillators constituting the crystal, can only 
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interact through the anharmonic term of the potential 
leading to three or more phonon collisions. These inter- 
actions preserve neither the number of phonons nor their 
frequency in the collision process. Nevertheless, these 
three or four phonons interactions tend to restore thermal 
equilibrium i.e. to help the phonons to follow an equi- 
librium distribution function which can be easily deter- 
mined from thermodynamic equilibrium considerations. 
Thus, much of the studies modelize the collision term in 
the BTE by the so-called Relaxation Time Approxima- 
tion : the distribution function /(r, p, i) relaxes to an 
equilibrium function /°(r, p) on a time scale r(p). The 
BTE resulting from this approximation is nothing but the 
Radiative Transfer Equation without scattering^. All 
the numerical tools developed in thermal radiation can 
therefore be used in this case. The key point in this model 
is to calculate a suitable r(p) in order to characterize the 
collisions. 

In the middle of the twentieth century, a great theo- 
retical effort has been done to modclizcd the relaxation 
times of the phonons in a bulk material. At ambient 
temperature, it has been shown that the main contri- 
bution to the relaxation time finds its origin in the an- 
harmonic phonon interaction. Two kinds of processes 
can be identified. The so-called Normal processes (iV) 
which maintain the momentum in the collision and the 
Umklapp processes (U) which do not preserve the mo- 
mentum. The former do not affect the material ther- 
mal resistance contrary to the latter. These Umklapp 
processes follow selection rulesi^ and it is an amazing 
feat to calculate themi^. Actually, the relaxation times 
are known for a small number of materials and often 
in the bulk situations. In the case of semi-conductors 
such as silicon (Si), germanium (Ge)i£*i& and gallium ar- 
senide (GaAs)i£ii&, the relaxation times have allowed to 
compute semi- analytically thermal conductivities in good 
agreement with measurements. Resolution of the BTE 
have been achieved on these materials in bulk situations, 
thin film or superlattice configurationi2ii2i22i2i. At short 
time scale, these resolutions have been compared to clas- 
sical solutions 2 ^ 2 ^ and some modifications of the BTE 
have been proposed-! The resolutions based on the Dis- 
crete Ordinates Method or on the Finite Volumes Method 
converge very quickly numerically but have a major draw- 
back : they are governed by a single relaxation time tak- 
ing into account all the different processes of relaxation 
such as the anharmonic interactions between phonons, 
the interactions with impurities and dislocations or the 
scattering on the material boundaries. The Matthiesen 
rule which stands that the inverse of the total relaxation 
term is the sum of the relaxation times due to every dif- 
ferent phenomena is usually used. In the context of the 
BTE in the relaxation time approximation, this means 
that all the different interaction or scattering phenom- 
ena tends to restore thermal equilibrium. 

An alternative way to solve the BTE is the Monte 
Carlo method. This method is quite computer time 
greedy because it necessitates to follow a large number 



of energy carriers, but it becomes competitive when the 
complexity of the problem increases, particularly for non- 
trivial geometries. This method is therefore useful in or- 
der to calculate the heat transfer in electronic devices of 
any shape. Moreover, in this method, different scattering 
phenomena (impurities scattering, boundary scattering 
and inelastic scattering) can be treated separately. The 
resolution of the BTE by the Monte-Carlo method has 
been performed for electrons 25 ' 26,27,28,29 '— but has been 
little used in the case of phonons. Peterson^ 1 performed 
a Monte Carlo simulation for phonons in the Debye ap- 
proximation with a single relaxation time. He presented 
results both in the transient regime and in equilibrium 
situation. Mazumder and Majumdar— followed Peter- 
son's approach but included in their simulation the dis- 
persion and the different acoustic polarization branches. 
They retrieved both the ballistic and the diffusion situa- 
tion but did not show any result in the transient regime. 
Another limit of this last paper is that the N processes 
and the U processes are not treated separately although 
they do not contribute in the same way to the conduc- 
tivity. 

The starting point in our work are these two contribu- 
tions. We follow individual phonons in a space divided 
into cells. The phonons after a drift phase are able to 
interact and to be scattered. The speed and the rate at 
which phonons scatter depends on the frequency. We en- 
sure that energy is conserved after each scattering pro- 
cess. This procedure is different whether the phonons 
interact through a N process or a U process. This paper 
is therefore an improvement of existing phonon Monte 
Carlo methods and is validated on simple examples such 
as a semiconductor film heated at two different temper- 
atures. 

Section II recalls the basic hypothesis governing the 
BTE. Fundamental quantities such as the number of 
phonons, the energy and the density of states are also de- 
fined. The phonons properties are also presented through 
their dispersion relations. Section III exposes the Monte 
Carlo method used in this paper. Boundary conditions, 
phonons drift and scattering procedures are given in de- 
tails. Section IV present transient results in the diffusion 
and ballistic regimes. Thermal conductivities of silicon 
and germanium between 100 K and 500 K are numerically 
estimated. Influence of conductivity thermal dependence 
on heat conduction within a slab is studied. Finally, size 
effects on phonons transport at very short scales are con- 
sidered. 



II. THEORY 

A. Boltzmann Transport Equation 

The Boltzmann Transport Equation (BTE) is used to 
model the phonons behavior in a crystal lattice. This 
equation is related to the variation of the distribution 
function f(t, r, K) which depends on time t, location r 
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and wave vector K. f(t, r, K) can also be defined as the 
mean particle number at time t in the <i 3 r volume around 
r with K wave vector and d 3 K accuracy. In the absence 
of external force, the BTE expression isS 



df _ df 



(1) 



collision 



with the phonon group velocity v g = Vkw. 

Integration of the distribution function over all the 
wave vectors of the first Brillouin zone and all the loca- 
tions leads to the phonons number N(t) at a given time 
in the crystal. The lhs term of Eq.QJ accounts for the 
phonons drift in the medium and the rhs term for the 
equilibrium restoration due to phonons collisions with 
themselves, impurities and boundaries. 

The collisional term modeling is the key point in the 
BTE resolution. In the case of photons, it can be quite 
easily modelized by an absorption term, an emission term 
and an elastic scattering term in which a scattering phase 
function relates a photon in the incoming and outgoing 
propagation directions during a scattering evenfcS. How- 
ever, in the case of phonons, there is no absorption, nor 
emission but only scattering events. Scattering events 
at the borders can simply be treated during the drift 
phase i.e. when a phonon reaches a border. Scattering 
with impurities can be treated similarly to the isotropic 
scattering of photons in thermal radiation. Scattering 
of phonons due to the anharmonic terms of the poten- 
tial are quite difficult to modelize. We know nevertheless 
that these terms are responsible of the thermal conduc- 
tivity i.e. tend to restore thermal equilibrium. Therefore, 
in this work we use the Relaxation Time Approximation 
for three phonons scattering processes. The collision time 
used in this formalism comes from Normal and Umklapp 
relaxation times which arc further estimated. 



E is the material volumic energy. It is obtain by sum- 
mation in Eq. J2J) of each quantum hw over the two polar- 
izations for transverse, longitudinal and optical modes of 
phonons propagation. Assuming that the phonons wave 
vectors are sufficiently dense in the K-space, the summa- 
tion over K can be replaced by an integral. Moreover, 
using D p (oj) the phonon density of state, we can achieve 
the integration in the frequency domain. This two mod- 
ifications yield 



- ) hwD p (uj)g p duj 
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with D p (uj)duj the number of vibrational modes in the 
frequency range [u>, uj + du>] for polarization p and g p the 
degeneracy of the considered branch. In the case of a 
three dimensional crystal (V — L 3 ) we have 



D p (uj)duj 
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The \ term in Eq.||2J) is the constant zero point energy 
which do not participate to the energy transfer in the 
material, therefore it has been suppressed. Using the 
group velocity definition, Eq.Q might be rewritten 
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The numerical scheme we are going to present is mainly 
based upon energy considerations. The previous expres- 
sion Eq.JSJ) will be also used to estimate the material 
temperature by means of a numerical inversion. 



B. Lattice modeling 



C. Dispersion curves 



As it has been exposed previously, the thermal 
behavior of the crystal can be considered from the 
phonons characteristics (location, velocity and polar- 
ization) within the medium. They might be obtained 
through the BTE solution since the distribution function 
can be easily related to the energy and therefore to the 
temperature. Using an integrated distribution function, 
one can express the total vibrational energy of the crystal 
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where (riK.p) is the local thermodynamic phonons popu- 
lation with polarization p and wave vector K described 
by the Bose-Einstein distribution function 
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Only a few studies on that topic take into account dis- 
persion. Indeed, frequency dependence makes calcula- 
tions longer accounting for velocity variation. However, 
realistic simulation of phonons propagation through the 
crystal must take into account interaction between the 
different branches. Here optical phonons are not consid- 
ered because of their low group velocity : they do not con- 
tribute significantly to the heat transfer. Consequently 
only transverse and longitudinal branches of silicon and 
germanium are presented here (Fig. Ill Cjl . We have made 
the common isotropic assumption for wave vector and 
consider the [001] direction in K-space. For silicon, we 
used data obtained from a quadratic fi^2, whereas ger- 
manium experimental curves 3 ^ have been fitted by cubic 
splines. Phonons group velocity has then been extracted 
from this data. Note that in silicon an germanium, two 
acoustic branches have been considered. The transverse 
branch is degenerate (gr = 2) whereas the longitudinal 
branch is non-degenerate (g^ = 1). 
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FIG. 1: Phonons dispersion curves for silicon and germanium 
in the first Brillouin zone, KmaxSr = 1.1326 • 10 10 m _1 and 
K maxGe = 1.1105- 10 10 m- 1 



III. MONTE CARLO METHOD 

The Monte Carlo technique has been widely used in 
order to solve transport equations. In the heat transfer 
field, Monte Carlo solutions of RTE are often considered 
as reference benchmarks. The method accuracy only lies 
on the number of sample used. Among others, the main 
advantages of this technique are : 

• The simple treatment of transient problems, 

• The ability to consider complex geometries, 

• The possibility to follow independently each scat- 
tering processes (for instance phonon-phonon, 
phonon-impurity and phonon-boundary processes). 

The main drawback is computational time. However 
this method remains a good choice between determinis- 
tic approaches such as the Discrete Ordinates Method 
(DOM) or "exact solution" such as those provided by 
molecular dynamic which is limited to very small struc- 
tures. 



A. Simulation domain and boundary conditions 

As it was said before, the geometry of the studied ma- 
terial does not matter. Here, a simple cubic cells stack 
(Fig. 1111 A(l is considered since it can be readily related 
to the plane wall geometry commonly used in thermal 
problems. Cylindrical cells or multidimensional stacks 
can be also considered in order to model nanowires or 
real semi-conductors. 

Concerning boundary conditions, we assume that the 
lateral walls of the cells (in x and y directions) are spec- 
ularly reflecting in most of the simulation cases. This 
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FIG. 2: Studied model. Phonons location, energy and veloc- 
ity are randomly chosen in each cell according to dispersion 
curves and local temperature. 



means that walls are adiabatic and perfectly smooth. 
Note also that in that case, the dimension in the x and y 
directions should not change the result in the simulation. 
Indeed, when reflection is specular on the lateral cells 
boundaries, the momentum is preserved in the z direc- 
tion. The heat flux and the temperature along z should 
thus not be affected. At both end of the medium, tem- 
perature is assumed to be constant. Therefore, energy in 
the first and the last cell is calculated from equilibrium 
distribution functions. Incoming phonons in these cells 
are thermalized at each time step. Consequently theses 
cells act as blackbodies. 

At this stage, an important point is the choice of the 
three discretizations : temporal, spatial and spectral. 
Spatial discretization is directly related to the material 
geometry : usually cells length are about L z ~ 100 nm 
for micrometric objects and can be smaller in the case 
of thin films or nanowires for instance. The time step 
choice depends on two parameters : the cell size and the 
group velocity at a given frequency. In order to consider 
all scattering events and to avoid ballistic jump over sev- 
eral cells, we state that the time step must be lower than 
At < L z /V g 

max • 

The spectral discretization is uniform, we used N 
1 000 spectral bins in the range [0,u)LAmax\- We have 
checked that larger discretizations do not increase the 
results accuracy. 



B. Initialization 

The first step of the simulation procedure, once 
medium, geometry and mesh have been chosen, is to ini- 
tialize the state of phonons within each cell describing the 
material. Hence, the number of phonons present in each 
cell is required. It will be obtained considering the local 
temperature within the cell and using a modified expres- 
sion of Eq. © . In this equation, energy is given for all 
the quanta fko associated to a spectral bin. Therefore, it 
can be rewritten to give the total number of phonons in 
a cell as 
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The number of phonons obtain with Eq. (JJJ is usually very 
large, for instance in a 10 nm silicon cube at 300 K, N can 
be estimated around 5.45 ■ 10 5 . In the case of nanoscale 
structures, direct simulations can be achieved if the tem- 
perature is relatively low. In the case of microscale sam- 
ples or multidimensional cell stacking, a weighting fac- 
tor shall be used to achieved Monte Carlo simulations. 
Hence, Peterson's^! technique has been used. The actual 
number of phonons N is divided by a constant weight W 
in order to obtain the number of simulated phonons N* 



W 



(8) 



In our simulations W maximal values are W ~ 10 4 for 
micrometric structures. 

During the initialization process, a temperature step 
is prescribed in the medium. The first cell being raised 
to the hot temperature Th, the last to the cold one T c . 
All the phonons in the intermediate boxes are also at 
T c . Associated theoretical energy in the whole structure 
is obtained from Eq.©. This energy should match the 
calculated energy E* within all the cells, written into the 
following form 



E* 



Ncell N* 

EE 



W x hw n 



(9) 



A new random number R is drawn: if R < Pla(^i) 
the phonon belongs to the LA branch otherwise it is a 
transverse one. 

The knowledge of the frequency and the polarization 
leads to the estimation of the phonon group velocity 
and the phonon wave vector merely using the dispersion 
curves and their derivatives. Assuming isotropy within 
the crystal, the direction f2 of the wave vector and the 
group velocity are obtained from two random numbers R 
and R' in the Cartesian coordinates (i, j, k) associated 
to lengths L x , L y and L z . Hence ft is written as 
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(13) 



The last operation of the initialization procedure is to 
give a random position to the phonon within the cell. In 
the grid previously considered, location of the n th phonon 
in the cell c is 



r c + L X R\ + L y R'j + L z R"k 



(14) 



where r c is the coordinates of the cell and R, R! and R" 
three random numbers. 



As a consequence, during the initialization, phonons 
should be added by packs of W at a given frequency, sam- 
pled from a normalized number density function F. Ac- 
cording to Mazumder and Majumdar work2£, this func- 
tion is constructed doing the cumulative summation of 
the number of phonons in the i th spectral bin over the 
total number of phonons Eq. Q) 



Fi(T) 



(10) 



In this process, a random number R is drawn (all the ran- 
dom numbers discussed here check < R < 1) and the 
corresponding value Fi gives the frequency u>i, knowing 
that Fi—\ < R < Fi location is achieved with bisection 
algorithm. The actual frequency of the phonon is ran- 
domly chosen in the spectral interval prescribing 



Vi = uj q ^ + (2R - 1) 



Alo 



(11) 



where wo, t is the central frequency of the i th interval. 

Once the frequency is known, the polarization of the 
phonon has to be determined. It can belong to the 
TA or LA branch with respect to the Bose-Einstein dis- 
tribution and the density of states. For a given fre- 
quency u>i, the number of phonon on each branch is: 
N LA (aji) = (n LA (u)i)) D LA (uji) and N TA (uji) = 2x 
{riT A (uJi)) DT A (oJi). The associated probability to find 
a LA phonon is expressed as 



N LA (LOi) + N TA {uj t 



(12) 



C. Drift 

Once the initialization stage is achieved, phonons are 
allowed to drift inside the nanostructure. Considering the 
time step At and their velocities, each phonon position 
is updated : r^rift = r i<j + v g At. In the case of shifting 
outside of the lateral boundaries (in i and j directions) 
the phonon is specularly reflected at the wall. In the case 
of diffuse reflection with a particular degree d (0 < d < 1 , 
d = purely specular, d — 1 purely diffuse) a random 
number R is drawn. When R is lower than d a new 
phonon propagation direction is calculated using Eq. 1|13[) . 

When a phonon reaches the bottom (z m i n ) or the top 
(zmax) of a cell, it is allowed to carry on its way in the 
previous or next cell respectively. As a result, it is go- 
ing to modify the cell energy and by extension its local 
temperature. At the end of the drift phase, the actual en- 
ergy E* is computed in all the cells using Eq.©. Then, 
the actual temperature T is obtained with Eq.© doing 
a Newton-Raphson inversion 36 . Phonons drifting in the 
first and last cellules are thermalized to the cold or hot 
temperature in order to keep boundary cells acting as 
blackbody sources. 



D. Scattering 

In the Monte Carlo simulation, the scattering pro- 
cess has been treated independently from the drift. The 
phonon-phonon scattering aims at restoring local ther- 
mal equilibrium in the crystal since it changes phonons 
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frequency. Collisions with impurities or crystal defects as 
well as boundary scattering do not change frequency but 
solely the direction fl. These last phenomena are signifi- 
cant when low temperatures are reached and the phonon 
mean free path becomes large. Here only three-phonon 
interactions have been considered. 

As already said before, there are two kinds of three- 
phonon processes : Normal processes (JV) which pre- 
serve momentum and Umklapp processes (U) which do 
not preserve momentum by a reciprocal lattice vector. 
These two mechanisms have consequences on the ther- 
mal conductivity of the crystal. When the temperature 
is sufficiently high (T > Toebye), U processes become 
significant and directly modify heat propagation due to 
the resistivity effect on energy transport. On the other 
hand, Normal scattering also affects heat transfer since 
it modifies frequency distribution of the phonons. For 
phonons described by (p, w, K) and (p',u/,K') scatter- 
ing to {p", u", K"), the following relations are checked 



energy : fko + tkj' <-> fko" 
N processes : K + K' <-> K" 
[/processes : K + K' <-> K" 



G 



(15) 



where G is a lattice reciprocal vector. Scattering also in- 
volves polarization in the way that acoustic transverse 
and longitudinal phonons can interact. According to 
Srivastava 3 ^ for N and U processes different combina- 
tions are possible 

N and U processes : T + T ^ L and T + L ^ L 
N processes only : T + T ^ T and L + L ^ L 

(16) 

For the N processes only, all the participating phonons 
must be collinear to achieve scattering. Usually these 
interactions are neglected. 

Direct simulation of phonons scattering is an awkward 
challenge. With Monte Carlo simulations, it is possible 
to modelize phonons collisions with neighbors as in the 
gas kinetic theory calculating a three particles interac- 
tion cross section. However, in the present study, the 
frequency discretization might not be sufficiently thin to 
assess every three phonons processes. Thus the collisional 
process is treated in the Relaxation Time Approximation. 
Several studies on that topic have been carried out since 
the early work of Klemensi 3 -, a detailed paper of Han and 
Klcmens 3 ^ recalled them. 

Relaxation times r have been proposed for several crys- 
tals. They depend on the scattering processes, the tem- 
perature and the frequency. Holland's work on silicon^ 
and the recent study of Singh for germanium 3 ^ provide 
various r values. The independence of the scattering 
processes is used to consider a global three phonons in- 
verse relaxation times accounting for N and U processes 
tnu- It nas been obtained using the Mathiessen rule 

( t nu = t n + T u _ ) ■ 

In order to be implemented in the Monte Carlo simu- 
lation, the scattering routine requires an associated col- 
lision probability P sca t- This one is derived saying that 



the probability for a phonon to be scattered between t 
and t + dt is dt/r. Thus, 



P.,, 



1 — exp 



-At 

tnu 



(17) 



A random number R is drawn, if R < P sca t the phonon is 
scattered. As a result new : frequency, polarization, wave 
vector, group velocity and direction have to be resample 
with respect of energy and momentum conservation. 

In previous studies on that topic; 31 ! 32 the frequency 
sampling after collision was achieved from the normal- 
ized number density function F at the actual tempera- 
ture T of the cell obtained at the end of the drift proce- 
dure. In this approach the actual energy after the scat- 
tering stage is usually different from the "target" one 
obtained with temperature T. Hence a subsequent "cre- 
ation/destruction" scheme is necessary to ensure energy 
conservation. In fact, in the preceding procedure at ther- 
mal equilibrium, the probability of destroying a phonon 
of frequency lo and polarization p is different from the 
probability of creating this phonon. This means that the 
Kirchhoff law (creation balances destruction) is not re- 
spected. In order to create phonons at the same rate they 
are destroyed at thermal equilibrium, the distribution 
function used to sample the frequencies of the phonons 
after scattering has to be modulated by the probability 
of scatterring. So we define a new distribution function 



F S cat (T) 



scat j 



ESl^i(r) X Pscat j 



(18) 



Taking into account the scattering probability in the 
distribution function F sca t ensures that a destructed 
phonon on both transverse and longitudinal branches can 
be resample with a not to weak energy as it can be seen 

(Fig. meg. 
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FIG. 3: Normalized number density function in silicon with 
and without P sca t correction 
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According to the described simulation procedure af- 
ter initialization step, phonons in cell c are described by 
[T c , F{T C ), N*(T C ), E*{T C )}. They are allowed to drift 
and the state of cell c before scattering is [T c , F(T C ), 
N *(T C ), E *(T C )]. Then three-phonons collisions occur 
and change energy by frequency resetting of the colliding 
phonons (using the distribution function F scat , leading to 
the final state [f D , F scat (f c ), N'*(T C ), E"*(f c )]. Hence 
energy can be express as : 

e"* = ]T htJ, + huj t (19) 

i=l i=l 
^F scat (T c ) ^F(T C ) 

Furthermore the number of colliding phonons can sim- 
ply be express as N'* cat {f c ) = P scat x N'*{T C ). One 
then sees, that, if we want to preserve energy at thermal 
equilibrium, the new normalized number density func- 
tion F SC at must take into account this collisional proba- 
bility. Energy conservation during Monte Carlo simula- 
tion might be noticed from frequency distribution (Fig. 
IIII Dl that matches the theoretical distribution give by 
Eq.Q. Concerning momentum conservation, the task is 
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FIG. 4: Frequency spectra at 300K and 500K for silicon 

harder to address since the Monte Carlo process consid- 
ers phonons one by one. Consequently triadic N or U 
interactions cannot be rigorously treated. In a first ap- 
proach, we propose the following procedure to take into 
account the fact that U processes contribute to the ther- 
mal resistance whereas the N processes do not. When 
the phonons scatter through a U process, their directions 
after scattering are randomly chosen as in the initializa- 
tion procedure. Therefore, these phonons are randomly 
scattered and contribute to the diffusion of heat. On the 
contrary, it is assume that scattering phonons experienc- 
ing a N process do not change their propagation direction 

n. 

Statistically, for a given temperature and frequency, 
the phonons are destroyed by scattering at the same rate 



they appear. A phonon which scatters has a great chance 
to be replaced in the computation by a phonon of a near 
frequency. Therefore, by this treatment, the N processes 
"approximately" preserve momentum. Nevertheless, a 
more accurate treatment should be done in order to re- 
spect exactly the momentum in the N processes. For a 
plan-parallel geometry, it seems possible to guarantee the 
momentum conservation in a single direction. 

In fact, the relaxation time estimation^- states that 
there is a frequency limit uJu m it for the transverse acous- 
tic branch. uJu m it actually corresponds to K = ~K max /2. 
Below this limit frequency, there is no U processes. On 
the other hand for w > wumit, N processes are no longer 
considered and the propagation direction must be resam- 
pled in the case of a collision. In what concerns the longi- 
tudinal acoustic branch, there is no limit frequency. Ac- 
cording to Holland 16 only N processes exist. However 
applying this assumption implies that momentum has to 
be conserved for each scattering event involving a LA 
phonon. This leads to thermal conductivity values higher 
than the theoretical ones for temperatures between 100 
K and 250 K. In order to ensure a more realistic mo- 
mentum conservation we set that half of the colliding 
phonons keep their original $7, the others (U processes) 
are directionally resampled. 



IV. RESULTS AND DISCUSSION 

Different kinds of simulations have been performed so 
as to check the computational method. Tests in both 
diffusion and ballistic regimes are carried out for silicon 
and germanium. Moreover, if small thermal gradients are 
considered, one can estimate the thermal conductivity k 
from the heat flux through the structure. This has been 
realized for Si and Ge between 100 K and 500 K. 

Knowing that the conductivity varies with tempera- 
ture according to a power law in the case of Si and Ge 
for T greater than 100 K, it is obvious that a large ther- 
mal gradient applied to our media should not bring a 
purely linear solution. Hence simulations in this specific 
case have been done. The model ability to correctly pre- 
dict steady state has been confirmed with comparison to 
analytical solution. 

Eventually, we studied size effects on thermal behavior 
of nanostructures. It appears that the ballistic regime 
can be retrieved at room temperature when the sample 
size is close to the nanometer scale. 



A. High temperature transient calculations 

Concerning high temperature transient calculations, 
the simulated case is described by the following parame- 
ters : 

• Hot and cold temperatures: Th = 310 K and T c = 
290 K, 
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• Medium geometry: stack of 40 cellules (L x = L y = 
5 • 10~ 7 m, L z = 5- 10~ 8 m), 

• Time step and spectral discretization: At = 5 ps 
and N b = 1000 bins, 



• Weighting factor: W 
8 • 10 4 for Ge. 



3.5 • 10 4 for Si and W = 



Both materials were tested. Germanium calculation re- 
sults are presented here (Fig. IIV A(l . In order to assess 
the Monte Carlo solution, transient theoretical compari- 
son exists in the case of the Fourier limit. Nevertheless, 
it requires that the thermal diffusivity a remains con- 
stant. In the chosen temperature range, according to 
the IOFFE database 4 ^, Ge thermal diffusivity is equal to 
a = 0.36-10~ 4 m 2 s -1 . 

The considered test case has been described in Ozisjk's 
book 41 - on heat conduction equation. Within the de- 
scribed structure heat transfer is along z axis and an- 
alytical solution for one-dimensional medium could be 
obtained from integral transform. Temperature distribu- 
tion in the slab is given by an infinite sum that requires 
enough terms in the case of short time calculation. How- 
ever a simpler analytical solution might be obtained with 
Laplace's transform 



T(z,t)-T(L,t) 
T(0,t)-T(L,t) 



erfc 
erfc 



2Vcrf 
2L + z 



erfc 



(2L-z 
\ 2y/od 



2Vat 



(20) 



with erfc the complementary error function. The theoret- 
ical solution is only valid for short time and its accuracy 
is better than 1% if the Fourier number (Fo = at/L 2 ) 
check Fo < 0.7. In the case of a 2/zm germanium slab 
it leads to t < 78 ns, which is largely enough to reach 
steady state. 
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FIG. 5: Transient temperature in Fourier's regime for germa- 
nium and comparison with the analytical solution of heat con- 
duction equation with a constant thermal diffusivity (ptGe = 
0.36 • 10~ 4 m 2 s~ 1 ) (dotted curves) 



The calculated values have been obtained from ten sim- 
ulations averaged Fig. IIV Al the random numbers seed 
being reset for each computations. Monte Carlo model 
ability to predict correctly temperature from the first mo- 
ments till steady state is clearly illustrated. The remain- 
ing noise can be reduced with lower values of W weighting 
factor. Diffusion regime is obtained after 30 ns. Similar 
results are obtained for silicon, however the modelized 
slab has to be larger (L — 4^m) because ballistic effects 
are observed near the cold limit. This point will be dis- 
cussed latter. 



B. Low temperature transient calculations 

For low temperatures, heat transport inside the slab is 
different since phonons interactions change. In this case 
U collisions are neglectibles and the only resistive pro- 
cesses are assigned to impurities, defects and boundary 
scattering. These phenomena have to be carefully exam- 
ined in the case of thermal conductivity estimation be- 
low 100 K. In fact, for very low temperatures the phonon 
mean free path grows and becomes larger than the struc- 
ture length. Hence, phonons can travel from hot to cold 
extremity without colliding. This is the ballistic regime 
similar to the one observed with photons exchanged be- 
tween two black plates at different temperatures 4 ^. In 
this peculiar case temperature in steady state is equal to 
the following constant value 



Tba. 



llistic 



T 4 + T 4 



1/4 



(21) 



The simulation case parameters are 



• Hot and cold temperatures: T/,. = 11.88 K, T c = 3 

K and T baUistic = 10 K, 

• Medium geometry: stack of 40 cellules (L x — L y — 
5 • 10~ 7 m, L z = 2.5 • 10" 7 m), 

• Time step and spectral discretization: At = 5 ps 
and N b = 1000 bins, 

• Weighting factor: W = 20 for Si and W = 30 for 
Ge. 



Results for silicon and germanium (FigflVB} give the 
expected results for the ballistic limit. It can be noticed 
that the current representation exhibits an artificial link 
between black boundaries and the first medium cell due 
to the spatial discretization. It can be seen that hot 
phonons do not fly straight toward the cold limit. More 
than 1 ns is necessary to heat the last cell in the case of 
silicon. This is in agreement with velocities prescribed 
by dispersion curves. Heat propagation in germanium is 
slower since phonons group speed is also lower. Results 
at low temperatures obtained with our method have been 
already be predicted by Joshi and Majumdar— in similar 
cases, who applied successfully the Equation of Phonon 
Radiative Transfer (EPRT) in ballistic regime. 
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C. Si and Ge thermal conductivities 

There are several ways to perform the thermal con- 
ductivity calculation of a semiconductor. Among these 
techniques, Holland's method— based on phonons kinetic 
theory was largely employed. Molecular dynamic simu- 
lations can be also used to obtain k. In the present study 
thermal conductivity has been determined knowing the 
heat flux (phonon energy transport) through the medium 
for a given thermal gradient directly applying Fourier's 
law. As in Mazumder's workS, the temperature differ- 
ence between hot ans cold extremities is set to 20 K so 
as to determine average conductivities. The phonon heat 
flux is calculated along z axis according to the following 
relation 



JY* 



= Whuj n \r g n • k 



(22) 



is large enough to assume the acoustic thick limit. Fur- 
thermore, this calcul benefits from recent relaxation time 
estimation which have been fixed with a good precision^. 
The influence of these factors on calculated conductivity 
is usually strong. Silicon results are also close to the bulk 
ones till 150 K where the relative error is equal to 7%. 
For lower temperatures, discrepancy between theory and 
simulation increases. This gap can be assigned to size 
effects since phonon mean free path grows when temper- 
ature is falling. Here, it becomes similar to the slab size. 
Yet, if we refer to Asheghi^ work on thin films thermal 
conductivity, at temperatures below 100 K, significantly 
decreases in comparison with bulk due to stronger reduc- 
tion of phonon mean free path by boundaries. Actually 
for pure 3 /iin silicon film thermal conductivity is close to 
600 Wm-'K- 1 at 100 K^i. This value is comparable to 
the 658 Wm-'K" 1 obtained for our 2 /jm film by Monte 
Carlo simulation. 



Simulations have been carried on between 100 K and 
500 K fFig llV Cjl on 2 fim thick samples. 

Comparison of the Monte Carlo calculated conductiv- 
ities is achieved with bulk data. Solid and dash curves 
are linear power law regression of theoretical data in the 
considered thermal range. These values are used in next 
analytical calculations. We have appraised for 200 K 
< T < 600 K 

exp (12.570) 
X s* i 1 ) - — ^326 — 

exp (10.659) 
X Ge (2 ) (23) 

For germanium a very good agreement is obtained with 
bulk values in the whole temperature domain. The max- 
imum relative error being under 8%. In this case, at 100 
K, phonon mean free path is lower to the micrometer 
according to Dames^. Consequently the stucture size 



D. Effect of non linear conductivity 

In this fourth part, transient simulations with samples 
heated under a large thermal step have been conducted. 
The purpose of such calculations was to underling the 
model capacity to correctly predict steady state when 
medium properties (fc(T)) vary with temperature. In the 
previous part thermal conductivities of both bulk materi- 
als have been estimated with a power law Eq. (|23|l . Hence, 
the analytic solution for temperature profile within a slab 
can be easily determined in steady state by the resolution 
of a first order differential equation as 



^ Steady state 



L 



(7 + 1) 



1 L) h 



1/(7+1) 



(24) 



where conductivity can be written as A(T) — C x T 7 . 
In order to avoid boundary effects in the case of silicon, 
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4/xm thick sample, with larger cells (L z = 1 • 10~ 7 m), is 
used. The initial geometry is kept for germanium slab. 
Temperatures are now : Th — 500 K and T c — 250 K. 
Time step remains equal to 5 ps. In both cases fFia llV D|) 
Monte Carlo simulations give a very good estimation of 
the steady state behavior. Results are averaged over five 
computations on the last 1000 time steps (i.e. after 45 ns 
of elapsed time) . At the cold limit of germanium sample a 
weak deviation exits between simulation and theory. The 
relative error on temperature remains smaller than 1.5%, 
in this area. This mismatch could be assigned to bound- 
ary effects, associating diffusive and ballistic regimes near 
the limits as we will detailled it in the last part. It could 
also be due to the accuracy of the bulk thermal conduc- 
tivity fitting at low temperatures. 




zftim) z(|im) 



FIG. 8: Steady state temperature in Fourier's regime for sil- 
icon and germanium in the case of a large thermal gradient; 
comparison to the heat conduction equation analytical solu- 
tion for temperature dependent conductivity. 

Besides, inversion of such curves can theoretically pro- 
vide the variation of A; on a given thermal range, as long 
as the medium is in the acoustic thick limit. 



Weighting parameters and lateral cells lengths are ad- 
justed in order to keep approximatetively 18000 phonons 
in each cell. 

Temperature profiles when steady state is reached have 
been plotted for each sample versus adimensionnal length 
z/L. Comparison to diffusive and ballistic regime is dis- 
played in Fie llV El With the imposed boundary temper- 
atures the ballistic limit is equal to Tb a m s ti C — 300.5 K. 
In the case of structures length lower than 200 nm ballis- 




0.2 0.4 0.6 0.8 1 

z/L 



FIG. 9: Steady state temperature for silicon, influence of the 
slab thickness; comparison to the analytical solution in the 
diffusive and ballistic limits. 

tic trend mixed to phonons diffuse transport is observed. 
The temperature profile gets closer to the ballistic limit 
for sample size around the nanometer scale. Neverthe- 
less, this approximately represents ten atom layers and 
therefore might encounter the modelization limitation. 
On the contrary, in a silicon sample thicker than 4 /xm, 
temperature reaches the Fourier's regime and can be sim- 
ilarly obtained with heat conduction equation at least 
cost. 



E. Size effect on heat diffusion 

From the previous calculations, it is obvious that the 
phonon mean free path modification with the tempera- 
ture acts as a major factor in heat conduction. So, if 
the structure size is adjusted in order to match the mean 
free path at any temperature, ballistic phenomena should 
be observed. In this study only silicon is used and the 
simulation parameters are 

• Hot and cold temperatures: = 310 K and T c = 
290 K, 

• Number of cells: 40, 

• Total length and time step (L, At): (2 nm ,5 • 10~ 3 
ps), (20 nm ,5-10- 2 ps), (200 nm ,5-10 _1 ps), (2^m 
,5 ps) and (4/mi ,5 ps). 



V. CONCLUSIONS 

An improved Monte Carlo scheme that allows transient 
heat transfer calculations at time and space nanoscales, 
on the basis of phonons transport, has been presented. 
This model accounts for phonons transitions between lon- 
gitudinal and transverse acoustic branches and can be 
simply applied to several semiconductors if their disper- 
sion relations are known. A particular attention has been 
paid to the energy and momentum conservation during 
collision process. 

Numerical results forecast have been assessed in differ- 
ent heat transfer modes. In slab configuration, a good 
agreement was found for both extreme phonon motion 
which are the diffusive and ballistic ones. Bulk thermal 
conductivities of silicon and germanium have been nu- 
merically retrieved with a maximal error lower than 8%. 
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Besides, our Monte Carlo model correctly predict tem- 
perature profile in more peculiar situations, when strong 
thermal gradient or very small sizes are encountered. 

Nevertheless some key points need to be refined. 
Among them momentum conservation procedure might 
be improved, especially for one dimensional applications. 
Optical phonons effect on heat transport were neglected. 
However according to the recent study of Narumanch: 23 
they must be taken into account especially for capaci- 
tive properties prediction. Regarding the collision pro- 
cess, improvements might be expected. Using theoret- 
ical values of r recalled by Han and KlemensS, direct 
calculation of phonon scattering relaxation time can be 
realized in each authorized spectral bin. Hence, a more 



realistic approach of three phonons interactions should 
be achieved. 

We are currently working on this improvements but 
also on other potential implementation of the method 
such as those related to the superlattices and the 
nanowires. 
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